On the complexity of solving initial value problems 
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ABSTRACT 

In this paper we prove that computing the solution of an 
initial-value problem y = p{y) with initial condition y{to) — 
yo £ R"^ at time to + T with precision e~'^ where p is a 
vector of polynomials can be done in time polynomial in 
the value of T, n and Y = supt^^^^^jn Contrary 
to existing results, our algorithm works for any vector of 
polynomials p over any bounded or unbounded domain and 
has a guaranteed complexity and precision. In particular 
we do not assume p to be fixed, nor the solution to lie in 
a compact domain, nor we assume that p has a Lipschitz 
constant. 

1. INTRODUCTION 

Solving initial-value problems (IVPs) defined with ordi- 
nary differential equations (ODEs) is of great interest, both 
in practice and in theory. Many algorithms have been de- 
vised to solve IVPs, but they usually only satisfy one of the 
following conditions: (i) they are guaranteed to find a solu- 
tion of the IVP with a given precision (accuracy) (ii) they 
are fast (efficiency). It is not easy to find an algorithm which 
satisfies both (i) and (ii) since, in practice, efficiency comes 
at the cost of accuracy and vice versa. 

Observe in particular that when dealing with functions de- 
fined on an unbounded domain (say R, the set of reals), nu- 
merical methods with guaranteed accuracy for solving ODE 
are not polynomiaQ even for very basic functions. Indeed, 

*In the classical sense, th at is to say to be formal, in the 
sense of recursive analysis WeiOO . 
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usual methods for numerical integrations (including basic 
Euler's method, Runge Kutta's methods, etc) fall in the 
general theory of n-order methods for some n. They do 
require a polynomial number of steps for general functions 
over a compact domain [a,b], but are not always polynomial 
over unbounded domains: computing f{t), hence solving the 
ODE over a domain of the form [0, t] is not done in a number 
of steps polynomial in t without further hypotheses on func- 
tion / by any n-order method (see e.g. the general theory in 
[Dem96] ). This has been already observed in [SmiOGl, and 
ODEs have been claimed to be solvable in polynomial time 



SmiOG for some classes of functions by using methods of 



order n with n depending on t, but without a full proof. 

Somehow, part of the problem is that solving y — f(y) 
in full generality requires some knowledge of /. Generally 
speaking, most algorithms only work over a fixed, specified 
compact domain, since in this case / is usually Lipschitz (it 
is well-known that any function over a compact is also 
Lipschitz there). 

In this paper we present an algorithm which has been 
designed to work over an unbounded domain, be it time 
or space. More precisely, we do not assume that the solu- 
tion lies in a compact domain. Our algorithm guarantees 
the precision (given as input) and its running time is ana- 
lyzed. Achieving both aims is especially complicated over 
unbounded domains, because we do not have access to the 
classic and fundamental assumption that / is Lipschitz. 

The proposed method is based on the idea of using a 
varying order: the order n is chosen accordingly to the re- 
quired precision and other parameters. Compared to [Smi06| 
(which takes a physicist's point of view, being more inter- 
ested on the physics of the n-body problem than on studying 
the problem from a numerical analysis or recursive analysis 
perspective), we provide full proofs, and we state precisely 
the required hypotheses. 

Even though many useful functions / are locally Lipschitz 
functions, there is an inherent chicken-and-egg problem in 
using this hypothesis. We can always compute some approx- 
imation z{t) of the correct solution y{t). However, to know 
the error we made computing z{T), one often needs a local 



Lipschitz constant L which is valid for some set A with the 
property that y{t),z{t) £ A for all t £ [to,T]. However, we 
can only obtain the Lipschitz constant if we know A, but we 
can only know A from z{t) if we know the error we made 
computing z{t), which was the problem we were trying to 
solve in the first place. 

In this paper we are interested in obtaining an algorithm 
which provides an approximation of the solution oiy = f{y), 
bounded by error e, where e is given as input (as well as other 
parameters). We will then analyze the impact of the various 
parameters on the runtime of the algorithm. 

We mainly restrict in this paper our analysis to the sim- 
ple yet broad class of differential equations of the form |[T|, 
that is y = p{y) where p is a vector of polynomials. This 
is motivated first by the fact that most ODEs using usual 
functions from Analysis can be rewritten equivalently into 
the format ([T]) 



GCB08 . We denote the solutions of 



this kind of problems as PI VP functions. 

The necessary input parameters for the algorithm will be 
specified as follows; Given p and a solution y : J — >■ R'' to 
the PIVP ([T]), we want to compute y(T) with a precision 
2"'" knowing that Vt T, \\y{T)\\^ < Y. Our parameters 
are thus d the dimension of the system, T the time at which 
to compute, jj, the precision, Y a bound on the solution and 
k — deg(p) the degree of the polynomials. 

Our second motivation for studying the particular class 
([T]) of polynomial IVPs comes from the study of the General 
Purpose Analog Computer (GPAC) [Sha41| . The GPAC is 
an analog model of computation introduced by Claude Shan- 
non as an idealization of an analog computer, the Differential 
Analyzer, which most well-known implementation was done 
in 1931 by Vannevar Bush [Bus31| . Differential Analyzers 
have been used intensively up to the 1950's as computa- 
tional machines to solve various problems from ballistic to 
aircraft design, before the era of digital computations that 
was boosted by the invention of the transistor |Wil96 . 

It is known that any GPAC can be described equivalently 
as the solution of a PIVP [GC03 . It has been shown that 



the GPAC (and thus PIVP functions) are equivalent to Tur- 
ing machines from a computability point of view |GCB08| , 
[BCGHOTa . However, it is unknown whether this equiva- 
lence holds at a complexity level. This work is a required and 
substantial step towards comparing the GPAC to Turing ma- 
chines, or if one prefers in more provocative terms, in proving 
the not clear fact that analog computation is not stronger 
than digital computation. See [BC081 and |BCGH07b| for 



more discussions. 

In other words, with the results presented in this paper, 
we seek not only to understand what is the computational 
complexity of the problem of solving PIVPs (which appear 
in many applications from Physics, etc.), but also to under- 
stand how analog and digital computational models can be 
related at a computational complexity level. 
Organization of the paper 

In Section [2] we introduce some notations and claim some 
basic results that will be useful later. In Section [3] we de- 
rive an explicit bound on the derivatives of the solution at 
any order. In Section |4] we derive an explicit bound on the 
divergence of two solutions given the initial difference. In 
Section |5] we apply the results of the previous section to de- 
rive an explicit error bound for a Taylor approximation of 
the solution at any point. In Section [6] we describe an algo- 
rithm to solve a PIVP and give an explicit complexity. We 



hence obtain the proof of our main result. 
Overview of the paper 

In order to compute the solution to a PIVP, we use a clas- 
sical multi-step method, but of varying order. At each step, 
we use a Taylor approximation of the solution at an arbi- 
trary order to approximate the solution. In Section [S] we 
explain how to compute the derivatives needed by the Tay- 
lor approximation. Furthermore, we need an explicit bound 
on the derivatives since our method is not of fixed order. 
Section [3] provides such a bound as a corollary. Since our 
method computes an approximation of the solution at each 
point, it will make slight errors that might amplify if not 
dealt with correctly. We can control the errors in two ways: 
by reducing the time step and by increasing the order of the 
method. A careful balance between those factors is needed. 
In Section |4] we explain how the error grows when the initial 
condition of the PIVP is perturbed. In Section[5]we quantify 
the overall error growth, by also including the error intro- 
duced by using Taylor approximations. Finally in Section [6] 
we put everything together and explain how to balance the 
different parameters to get our main result. 

2. NOTATION AND BASIC FACTS 

We will use the following notations: 



IK 



Xl, 



max \Xi\ 

l<i<n 



(2;i, . . . ,2;„)|| = Vl^t-l'"^ H + l^^np 

a; < y <4> Vi, Si < j/i 



id if n = 

j[n-i] Q J otherwise 



Safit)^ sup 11/(^)11. 



— /'"'(a) 
fc! 



it - af 



Note that is not an order on but just a notation. We 
will denote by W the principal branch of the Lambert W 
function which satisfies x = W{x)e^^^^ and W —1. 
We will consider the following ODE: 



(1) 



y =p{y) 
y{to)= yo 

where p : R'' — 5> R'' is a vector of polynomial. If p : R"* 
is polynomial, we write: 

p{Xi,...,Xd) = 



where k is degree of pi that will be written deg(pi); We write 
|a| = Qfi + ■ • ■ + ad- We will also write: 

If p : R'' — )■ R'* is a vector of polynomial, we write deg(p) = 
max(deg(pi), . . . , deg(pd)) and Ep = max(Epi, . . . , Epd). With 
the previous notation, the following lemmas are obvious. 



Lemma 1. For any Q : R" — ^ R and any x G R", 
|Q(a;)| ^EQmax(l,||a;||l^^('3)) 

Lemma 2. For any polynomial Q ; R" — > R, a G N'' and 
X G R", if \a\ ^ deg(Q) then 

|g("'(x)i < ia|!EQmax(l,||a;|ll^^(«-l"l) 

3. NTH DERIVATIVE OF Y 

Given the relationship between y and p(j/) it is natural 
to try to extend it to compute the nth derivative of y. Of 
particular interest is the following remark: the derivatives 
at a point t only depend on y(t). Since an exact formula is 
difficult to obtain, we only give a recursive formula and try 
to bound the coefficients of the relationship obtained. 

Notice that since pi is a polynomial, it is infinitely dif- 
ferentiable and thus the partial derivatives commute. This 
means that for any a G N"^, p^"' is well-defined and is inde- 
pendent of the order in which the derivatives are taken. 

Proposition 1. Ify satisfies ([TJ for any t G I anddeg{p) 
k, define: 

r = io,fcf A = ii,d]xr 



V{t)={p^-\y{t)))^^ 



t G / 



Q)eA 



Then 



yt G /,Vn G ^(t) = Q^AV{t)) 

where Qi^n is a polynomial of degree at most n. Furthermore, 

Proof. First notice that Qi,„ has variables p^"\y{t)) so 
we will use dj^aQi,n to designate its partial derivatives. We 
will prove this result by induction on n. The case of n = 1 
is trivial: 

y'^it) ^ My{t)) Q,,i(V(t))=P.(y(i)) 

deg(Q,,i) = l = 1 = 0!''') 

Now fix n ^ 1. We will need /ij, G N'' such that {fik)i ~ Si^k- 
Elementary differential calculus gives: 



dy, 



(n) 



d 



dt dt 
d 



= EE {J2y''ii)9kP^"\y(m a.,«o.n(i^(t)) 



j=i Qer \fc=i / 
= EE (Ep^(2/W)p^"^"''(2/W)) d,.^QUV{t)) 

j = l agr \fc = l / 

Since Qi,„ is a polynomial, this proves that Qi,n+i is a poly- 
nomial. Furthermore a close look at the expression above 
makes it clear that each monomial has degree at most n + 1 
since every monomial of dj,aQi,niV{t)) has degree n—1 and 
is multiplied by the product of two variables of degree 1. 



Now, we can bound the sum of the coefficients. We will 
first need to give an explicit expression to (5i,„ so we write: 



Recall that the variables of Qi^„ are p^"'' so /? G 
and Pi^a makes perfect sense. Then: 

j=i aer \fc=i / 

j=iaer \|/3Kn / 

d 

= EE'^ E 

d 

= rf E l^/slEE''^-" 

d 

= rf E MJ2\p\ 

= duTiQi^ri 
dn(n- l)!''*' 



□ 

Corollary 1. Ify satisfies ([T]) fort G / anrf deg(p) — k, 
then 

||y<"'(t)|| n!*'*' max (^1, fclEpmax (l, 

4. DEPENDENCY IN THE INITIAL CON- 
DITION 

When using multi-steps methods to solve ODEs or simply 
when doing approximation, we might end up solving an ODE 
like ([T| with a wrong initial condition. This will of course 
affect the result of the computation since even if we could 
compute the solution with an infinite precision, the results 
would be different because of the dependency in the initial 
condition. For this reason we would like to evaluate this 
dependency numerically. Assuming y satisfies Q, we define 
the functional "I> as: 

$(io,yo,t) = y{t) 

Notice that the dependency in j/o is implicit but that in 
particular, $(to,yo,io) = yo- Also notice that $ and yo are 
vectors so we'll study the dependency of $i in yoj. 
We first recall the well-known Gronwall's inequality. 

Proposition 2 (Generalized Gronwall's inequality). 
Suppose ijj satisfies 

ip{t)^a{t)+[ f3{s)i;{s)ds, te[0,T] 
Jo 

with a(t) G R and /3(s) ^ 0. Then Vt G [0,r], 

V-W < a(t) + J a(s)/3(s)exp I3{u)du^ ds 



//, in addition, a is a non-decreasing function on [0,T], then 
ipit) a{t) exp /?(s)ds^ , t e [0, T] 

In order to apply Gronwall's inequality to (j), we will need 
to bound the Lipschitz constant for a multivariate poly- 
nomial. So consider a polynomial P G R[Xi, . . . ,Xd] and 
write: 

We first prove a lemma on monomials and then extend it 
to polynomials. 

Lemma 3. If a,b e R"^, a G N'* and ||a||^ , ||b||_^ sC M 
then: 

16" - a"| |q|MI°I"' lib- all 

I I \ I I II II oo 

Proof. One can see by induction that: 

i — l j <i j>i 

Since it is well know that for any integer n: 

n-l 

— a = (b — a) y a b 



i=0 



Thus we can deduce that: 

d 



i — 1 j-Ci j>i 
d ai-1 



3=0 



^||6-a|L^Ml"l-^a, 



«:|Q|||6-a|LM 



|q|-1 



□ 



We can use this result to obtain an explicit Lipschitz 
bound for the polynomial P: 

Lemma 4. For all a,b eR"^ such that \\a\\^ , \\b\\^ ^ M, 
\P{b) - P{a)\ ^ kM^~^T.P \\b - a\\^ 
where k = deg P. 
Proof. 

\Pib)-P{a)\ < Yl 

< |a„||Q|A/l"l-M|6-a|L 

s; kM"-^ ll&-a|L 



s; fcM'="^EP||6-a|| 



In order to evaluate the divergence between two solutions, 
we will need to solve a highly nonlinear equation (the so- 
called chicken-and-egg problem of the introduction). This 
lemma gives an explicit solution as well as an inequality 
result. 

Lemma 5. Let a,l3,x ^ and fc G N* then: 

X = ae^ <^ X = ae ^ ^ ' (2) 



^ has a solution and x ^ 4a 



Furthermore, 

k/3a' ^ i 

Proof. If a = this is trivial (a:: = 0), so we can assume 
that Of > and write a = e" and then a; > so we can write 
a; = e^. Then 

X — ae <^ X = a -\- pe 

<^ g = j3e^°'e^^ where g = x — a 

kg = Pe^^ where /3 = kpe^°' = fc/3a'^ 

<^ h = /3e'^ where h = kg 

~W{-P) 
<^ k{x - a) = -W{-kPa'') 

^x = a- yW(-kPa'') 
k 



<^ X = ae 



-^W{-k|3a'') 



□ 



And finally we can apply this result to <1>. 

Proposition 3. Let I = [a,b] and yo,zo e M''. Assume 
that y — "l>(o, 1/0,-) CLnd z = $(a, zq,-) '"^e defined over I. 
Let Y — SaV- Assume that \/t G /, 

||yo-zolLexp(fc^2'=Ep|t-a|(f + ^ i (3) 

Then \/t e I , 

II /j-MI ^ II II fc(2 + y(t))'=^lEp|t-a| 

where k = deg(p). 

Proof. Define M = max(y, Saz) and consider ij){t) = 
\\z(t) — y{t)\\^. By definition of $ we have: 

y{t) = yo+ p{y{u))du, t G / 

J a 

z(t) = zo + / p(z{u))du, t £ I 

J a 

Applying a few inequalities and Lemma |4] we get 

II 20 -yo 11^ + / \\p{z{u)) -p{y{u))\\^du 

J a 

= INo - yoll^ + / kM(t)'"''T,p'tP{u)du 

V ' Ja ^ ^ ' 



□ 



Finally, apply Proposition[2]witli a and /3 being non-decreasing 
functions to get 



ip{t) s; exp (^kMitf-'^Ep] 



t — a\ 



tel 



Then Vt G /, 



where k = deg(p). 



This inequality looks good except that there is an hidden 
dependency in M: M depends on ${a,yo,-) and <E>(a,j/o,-) 
and we seek one on $(a, j/o, •) only. Since tp is the difference 
between the two solutions, we have the following bound on 
M: 



M{t) s; Y{t) + Sai'i.f)^, tel (4) 
=*(t) 

Thus Vt G I, 

tPit) s; V(a) exp (^k {Y(t) + Ep|t - a|) 

G(y{t),*(t)) 

And since Y , ^l/ and G are non-decreasing functions, G{Y{t), '^(t)) 
is a non-decreasing function so we have: 

*(t)^G(y(t),*(t)), tei 

Consider the solution / : J — >■ M, J C j[^to: 
f{t) = ^(a) exp (fc2'= (y(t)'=-^ + fit)"^-') m - a!) (5) 

H(Y(t),f(t)) 

Since ([5| implies /(a) = = 'l'(a) and 5S x ^ y ^ 

Vz ^ 0,G{z,x) ^ H{z,y), then Vt G J, 'l'(f) sC /(f) and we 
can find an explicit expression for /: 



5. TAYLOR APPROXIMATION 

We first recall a simplified form of Taylor-Lagrange theo- 
rem which will be useful for our approximation step of the 
solution. 

Proposition 4 (Taylor-Lagrange). Leta,x eR, f e 
C''+^{[a,x]), then 



fix)-T^fix) < 



(^ - a)'' g r(k) 



k\ 



Sar>{x) 



fsf ^ f{t) = jp{a) exp (^fc2'=F(t)'=~'Sp| 



t — a\ 



X exp (fc2*Sp|t-a| /(*)'="') 



/3(t) 



^/(f) = a(t)exp (/3(f)/(t)^-^) 



Applying Lemmajsjand since ^ implies that (fc— l)/3(i)a(f)* ^ 
I we have: 

fit)^4:a(t), tel 

Notice that the case of fc = 1 is not handled by Lemma [5] 
but trivially gives the same result. Thus Vt G /, 

M{t) ^ Y{t) + 4iP{a) exp (^k2''Y{t)''~^Ep\t - aj^ 
^ 2 + Y(t) 

So finally, 

ip{t) sC TP{a) exp (^fe(2 y (t))'="^Epjf - aj^ , i G / 
□ 

Corollary 2. Let I = [a, 6] and yo,zo G R'*. Assume 
i/iffli y = $(o, 2/0,-) '"ic' -2: — $(a, zo,-) '"^e defined over I. 
Let Y = S'ay and M ^ | "^'^'^ assume that Wt G /, 

||yo-2o|Lexp('fc4'=Ep|t-al(l + y(f)'=-')) ^ (6) 



The idea is now to apply this result to the solution of Q 
and extend it in two directions: 

• Since the solution satisfies ([T]), we can use Corollary [l] 
to estimate the high-order derivative and the error 
bound. 

• Since our algorithm will make slight errors, we do not 
assume that we have the right initial condition; we 
want a general result with a perturbed solution and 
relate it to the expected solution using Corollary [2] 

Proposition 5. Let I = [a,b] andyo,zo G R''. Assume 
that y — ^{a,yo,-) and z = $(a, 2o,-) '^'^e defined over L. 
Let Y — Say and M 5; | and assume that Vt G /, 

|lj/o-2olLexp(fc4''-Ep|f-a|(l + y(f)'='^)) ^ (7) 

Then Vt G /, 

\\y{t) - T:z(t)\\^ ^ ^^+[d{t - a) (l + fe!Ep(l + + y(f))"))" 

where k = deg(p). 

Proof. Let A = \\y{t) - T^z{t)\\^, then 
A^\\yit)-zm^ + \\zit)-T:z{t)\\^ 

Apply Corollary |2] and Proposition [4] 



^ ^J.+ 



Saz'-^'it) 



Apply Corollary [T] 



^ ^J.+ 



max (1, fc!Epmax(l, Saz{t)'')y 



Use \\z{t)\\^ ^ lly(*)IL +M and apply Corollary|2] 



^ ^J.+ 



^ ^J.+ 



it - ")^ 

nl 

[t - af 



W}"> max (1, fc!Epmax(l, (fi + y(t)) ))' 
d!'") max (1, fc!Ep(l + + Y{t))'')" 



H + {d{t -a){l + fc!Ep(l + + y(f))'=))" 



^Notice that J can't be empty because a e J since /(a) 
4>{a)- We will see that J = L 



□ 



6. OUR MAIN RESULT 

First we need a lemma that will be helpful to compute the 
forward error. 

Lemma 6. Let a > 1 and b ^ 0, assume u € satisfies: 
u„+i ^ au„ +b, n ^ 

Then 

, a" - 1 

u„ ^ a uo + b — , n > 

a — 1 

Proof. By induction, the case n = is trivial and the 
induction step works as follows: 

Un+l < aUn + b 

^ n + 1 , I ^" ~ 1 t J, 

a ito + ab h 

a — 1 

^ n+i , , a(a"-l) + (a-l) 

^ a ^ uo + 

a — 1 

□ 



Proof. 



SC e 



E 



-£-A A 

< e e 



□ 



We now get our main result in a technical form: 

Theorem 1. //y satisfies ^ for t e I ^ [to,to + T], let 
k = deg(p), ne'M, T eQ+,Y eQ such that 

Then Algorithm^ above guarantees 

\\y{to + T) - SolvePIVP(to, yo, P, T,iu,N,iu)\\^^ e"^ 
with the following parameters 

M = (2 + Y)'= A ^ d{l + k\T.pM) iV = [TeA] 



T 

N 



B = fcrSpAM Lj^2 + fi + In(iV) + iVB 



Algorithm 1: NthDeriv 



input : The polynomial p of the PIVP 
input : The value z £ Q** of the function 
input : The order n of the derivative 
input : The precision ^ requested 
output: a: G Q'* 
1 Compute X such that \\x — i/'-"^(0) ^ where 

011 1 1 oo 



Algorithm 2: SolvePIVP 



input 
input 
input 
input 
input 
input 
output 
1 begin 

3 X -(r- yo 

4 for n 



The initial condition (to,2/o) G 
The polynomial p of the PIVP 
The total time step T G Q 
The precision ^ requested 
The number of steps A'' 
The order of the method uj 
a; G Q'* 



1 to iV do 

^bj-l A' 



E"=o' f^NthDerivCp, to + nA, x,iu,^ + A) 



Since at each step of the algorithm we compute an ap- 
proximation of the derivatives, we need a technical result to 
compute the total error made. We want to relate the er- 
ror between the value computed by the algorithm and the 
Taylor approximation, to the error made by computing the 
derivatives. 

Lemma 7. Let n,^,d G N*, A G Q+, 2,5 G (R'*)", as- 
sume that \\zi — ; 



SC e 



-e-A 



then 



1-1 



E 



fc! 



■2fc 



E 



fc! 



■Zk 



II - II ^ -JVS-ii-l 

Proof. Denote by t^ — to + nA and a;„ the value of x at 
the n^^ step of the algorithm. That is: 

y'"' = \\xr,+i ~ Ttlv^^'^tn + A)\\ «:e-" 

II 1 1 oo 

Notice that Xn+i is only an approximation of the Taylor 
approximation. We request an approximation up to e~'^, 
as a parameter of the algorithm. Lemma [7] ensures that 
this bound is indeed reached by computing the derivatives 
sufficiently precisely. Then define 

£n = \\Xn -yitn)\\^ 

By the choice of xo = yo, we have eo = l|yo — J/o|loo ^ 
g-jvB-fi-i ]\Jq^ assume that £„e'^ 5C |. After a few sim- 
plifications, Proposition [5] gives: 

En + l < £ne^ + (AA)" + e~" 

Now apply Lemma [6] 

p"^ - 1 



nB 

B^- 1 



Notice that 



n—l 



kB , nB 

e ^ ne 



Thus for n — N we get 

£jv ^ e e ^ -I- ((A^) +e j Aie 
By the choice of A we have 

AA ^ e"^ 
Thus by the choice of ui we have 



s; e 



+ 2e- 



-2-ii-ln(N)-NB jy^NB 



+ e 



Notice that we need to check that we indeed get a value 
smaller that i at the end otherwise we couldn't have applied 
Proposition [5] □ 



Lemma 8. // the coefficients of the vector of polynomial 
p are polynomial time computable, then for all z € Q"^ and 
n,^ G N*, NthDeriv(p, z, n, ^) has running time polynomial 
in the value of n and ^. 

Proof. From Proposition [T] we have: 

where Qi^„ is a polynomial of degree at most n and 
(n — From the proof, it is easy to see that the Qi,n 

are computable by induction in polynomial time in n and d 
since Qi,n is of degree n (thus has at most terms) and has 
the sum of its coefficients not larger that (n — (thus 
taking a space and time at most polylogarithmic in n to 
manipulate). Finally, V{t) can be computed with precision 
e~^ in time polynomial in ^ from the assumption on the 
coefficients of p. □ 

Corollary 3. // the coefficients of the vector of poly- 
nomial p are polynomial time computable, then for all to G 
Q, 2/0 G T e Q, AT, a; G N, SolvePIVP(to, yo,P, T, iV, uj) 
has running time polynomial in the value of T/N, ^, N and 

UJ. 

In less technical form: 

Theorem 2. There exists an algorithm A such that for 
any p vector of polynomial with polynomial time computable 
coefficients, yo G K'' polynomial time computable vector, 
to G Q, G N, T G Q and Y €Q such that Y > StoVito+T), 

\\A{p, yo, to, 11, T, Y) - ^{to,yo,to + r)||^ sC e"" 

Furthermore A{p, yo, to, fi, T, Y) is computed in time polyno- 
mial in the value of fi,T and Y. 

7. EXTENSION 

Our main result has the nice property that it requires 
minimal hypothesis and has maximal precision in its state- 
ment: all the parameters of the polynomials are kept. On 
the other hand, only considering differential equations of the 
form y — p{y) has two drawbacks: 

• it is not always possible; 

• even when possible, it might considerably increase the 
size of the system, and thus the complexity (which can 
be exponential is the size of the system) compared to 
a "non-expanded form" (e.g y = sin''''(j/) will expand 
to a 2fc-dimensional system). 

For this reason, we want to extend our result to differential 
equations of the form y — f{y) with restrictions on /. In- 
tuitively, our hypothesis will be the lemmas we had for p in 
the previous sections. That is, / will need to be computable 
quickly and its derivatives must not grow too fast. We only 
give the lemmas and proposition which are similar to the 
previous sections. The proofs are similar and we make a few 
remarks when relevant. For simplicity, we write the hypoth- 
esis only once, after introducing a small notation. Notice 
that in this section we assume d is a constant (since every- 
thing will be exponential in d anyway, we do not need to 
consider it). Hypothesis [Tlis about the differential equation 
verified by y. Hypothesis|2] is about value the derivatives of 
/. Hypothesis [3] is about the continuity modulus of / (or 



its "Lipschitz" constant). Finally, Hypothesis [4] is about the 
complexity of computing / and its derivatives. 

Notation: poly* {X #Y) = poly(X) 



poly(y) 



Hypothesis 1. Let I c R, J c R'' , y : I ^ J, f : J ^ 
to G Q, yo e K"^- Assume that 

y = f{y) 
yito)= yo 

Hypothesis 2. Assume that for all a G N'', x £ J , 
||/("'(x)j| ^poly* (||a;||^,|a|#|Q|) 

II I I OC 

Hypothesis 3. Assume that for all a,b e J, 

ll/W - /Wlloo ^ ll« - &II00 0(max||a||^ , ||6||^) 

where Q is a polynomial of degree k. 

Hypothesis 4. Assume that yo is polynomial time com- 
putable vector of reals and that f and its derivative are poly- 
nomial computable, that is, for all a & N, x € J, /'"'(x) 
is computable with precision e~^ in time polynomial in the 
value of\a\, ^ and ||2;||^. 

Hypothesis [l] will provides us a with a recursive formula 
for the derivatives of y, similarly to Proposition [l] Armed 
with Hypothesis |2] which more or less replaces Lemma [2] 
we can derive a bound on the derivative of y, similarly to 
Corollary [l] 

Proposition 6. Define 

Pn = {a G N"* I \a\ n} A„ = Jl, d] x r„ 



K(t) = (/i"'(y(t)))^^ 



tel 



(i,c<)6A„ 

Then 

Vi G J,Vn G N*,yf ^(t) = Q,,„(K-i(t)) 
where Qi^n is a polynomial of degree at most n. Furthermore, 
SQ,,„ (n- l)!'""' 

Corollary 4. 

\W"\t)\\ ^ poly* (lly(t)ll^, n#n) 

Similarly to Proposition [3j Hypothesis [3] will allow us to 
bound the divergence of two solutions given the initial dif- 
ference. We reuse the notation y — ^{a,yo,t) for / with its 
obvious meaning. 

Proposition 7. Let K — [a,b] and yo,zo G R''. Assume 
that y = $(a, t/o,-) and z = $(a, zo,-) are defined over K. 
Let Y — Say- Assume that Vi G /, 

\\yo - zoW^ S0(1 + Yit))" exp (EQ\t - a\{l + Yit))") ^ ^ 
Then Mt e I , 

The other results are basically the same except for the 
exact constant choices in the theorem. We get the same 
result at the end, that is y{t) is computable is polynomial 
time with respect to the same parameters. 
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